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Separated  Component-Based  Restoration 
of  Speckled  SAR  Images 

Vishal  M.  Patel,  Member,  IEEE ,  Glenn  R.  Easley,  Member,  IEEE ,  Rama  Chellappa,  Fellow,  IEEE , 

and  Nasser  M.  Nasrabadi,  Fellow,  IEEE 


Abstract — Many  coherent  imaging  modalities  such  as  synthetic 
aperture  radar  suffer  from  a  multiplicative  noise,  commonly 
referred  to  as  speckle,  which  often  makes  the  interpretation 
of  data  difficult.  An  effective  strategy  for  speckle  reduction  is 
to  use  a  dictionary  that  can  sparsely  represent  the  features  in 
the  speckled  image.  However,  such  approaches  fail  to  capture 
important  salient  features  such  as  texture.  In  this  paper,  we 
present  a  speckle  reduction  algorithm  that  handles  this  issue 
by  formulating  the  restoration  problem  so  that  the  structure  and 
texture  components  can  be  separately  estimated  with  different 
dictionaries.  To  solve  this  formulation,  an  iterative  algorithm 
based  on  surrogate  functionals  is  proposed.  Experiments  indicate 
the  proposed  method  performs  favorably  compared  to  state-of- 
the-art  speckle  reduction  methods. 

Index  Terms — Image  restoration,  multiplicative  noise,  speckle, 
synthetic  aperture  radar. 


I.  Introduction 

COHERENT  imaging  systems  such  as  synthetic  aper¬ 
ture  radar  (SAR),  holography,  ultrasound,  and  synthetic 
aperture  sonar  suffer  from  a  multiplicative  noise  known 
as  speckle  [1].  Speckle  appears  when  objects  illuminated 
by  coherent  radiation  have  surface  features  that  are  rough 
compared  with  the  illuminating  wavelength.  It  is  caused  by 
the  constructive  and  destructive  interference  of  the  coherent 
returns  scattered  by  many  elementary  reflectors  within  the 
resolution  cell.  Speckle  can  make  the  detection  and  interpre¬ 
tation  difficult  for  automated  as  well  as  human  observers. 
In  some  cases,  it  may  be  important  to  remove  speckle  to 
improve  applications  such  as  compression,  target  recognition, 
and  segmentation. 

Many  algorithms  have  been  developed  to  suppress  speckle 
noise  [2]-[7].  One  of  the  simplest  approaches  for  speckle 
noise  reduction  is  known  as  multilook  processing.  It  involves 
noncoherently  summing  the  independent  images  formed  from 
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L  independent  pieces  of  the  phase  history.  The  averaging 
process  reduces  the  noise  variance  by  a  factor  of  L.  However, 
this  often  results  in  the  reduction  of  the  spatial  resolution. 
Other  types  of  speckle  reduction  methods  are  based  on  spatial 
local  filtering  performed  after  the  formation  of  the  SAR  image. 
Various  filters  have  been  developed  that  avoid  the  loss  in  spa¬ 
tial  resolution  [2],  [8]— [1 1].  Some  of  these  methods  are  based 
on  a  window  processing  of  the  noisy  image.  Consequently, 
their  performance  depends  significantly  on  the  type,  direction, 
and  the  size  of  the  filter  used.  Furthermore,  some  of  these 
filters  often  fail  to  preserve  sharp  features  such  as  edges. 

To  overcome  some  of  these  limitations,  wavelet-based  meth¬ 
ods  are  often  utilized  [12]— [17],  in  which  noise  shrinkage 
is  applied  to  the  detailed  wavelet  coefficients  of  the  noisy 
image  [18],  [19].  Since  speckle  is  multiplicative  in  nature, 
some  of  these  methods  often  apply  the  logarithm  transform  to 
SAR  images  to  convert  the  multiplicative  noise  into  additive 
noise.  After  applying  soft  or  hard  thresholding  to  the  wavelet 
coefficients  of  the  logarithmically  transformed  image,  an  expo¬ 
nential  operation  is  employed  to  convert  the  logarithmically 
transformed  image  back  to  the  original  multiplicative  format. 

It  is  well  known  that  shrinkage-based  denoising  algorithms 
rely  on  the  sparsity  of  the  representation.  A  fixed  transform 
such  as  a  wavelet  transform  can  represent  a  piecewise  smooth 
image  sparsely  but  it  may  also  fail  to  represent  an  image 
with  textures  sparsely.  As  a  result,  the  overall  denoising 
performance  of  a  fixed  transform  on  an  image  containing  both 
piecewise  smooth  and  texture  components  can  be  inadequate. 

Another  popular  approach  for  restoring  speckled  images 
involves  total  variation  (TV)  regularization  [20],  where  the 
underlying  reflectivity  image  is  assumed  to  be  piecewise 
smooth  [21]-[23].  It  has  been  shown  that  TV  regularization 
often  yields  images  with  the  staircasing  effect  [24].  As  a  result, 
the  estimated  image  usually  contains  constant  regions,  and 
fine  details  such  as  textures  are  often  removed.  To  deal  with 
some  of  these  effects,  a  hybrid  method  that  uses  coefficient 
thresholding  and  TV  regularization  on  the  logarithm  of  the 
magnitude  image  or  the  log-image  domain  was  recently  pro¬ 
posed  [25].  In  particular,  hard  thresholding  on  the  curvelet 
coefficients  of  the  log-image  is  first  applied.  Then  a  variational 
method  that  uses  an  I  \  fidelity  to  the  thresholded  coefficients 
and  a  TV  regularization  in  the  log-image  domain  is  applied. 

Several  methods  have  been  proposed  that  use  a  combined 
dictionary  approach  to  image  restoration.  Suppose  that  we  are 
given  M  different  dictionaries  Dm,  m  =  1, ...  ,M;  then  one 
can  obtain  M  different  estimates  of  x  by  applying  either  hard 
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or  soft  threshold  to  the  coefficients  from  each  corresponding 
dictionary.  Let  x;  be  the  resulting  estimate  from  the  i  th 
dictionary.  Then,  a  simple  estimator  of  x  is  given  by  averaging 
M  individual  estimates 


x 


1 

M 


This  is  essentially  the  idea  behind  translation-invariant  wavelet 
denoising  by  cycle  spinning  [26],  where  the  final  estimate 
is  obtained  by  averaging  the  estimates  from  the  thresholded 
orthogonal  wavelet  coefficients  of  translated  versions  of  the 
original  image.  This  method  has  been  applied  for  speckle 
reduction  in  SAR  imagery  [12].  This  simple  method  suffers 
from  some  issues  in  practice,  as  it  weighs  equally  both 
good  and  bad  quality  estimates.  To  deal  with  this  problem, 
a  Bayesian  framework  to  optimally  combine  the  individual 
estimators  was  proposed  in  [27].  This  method  weights  each 
estimate  x;  at  each  sample  according  to  the  significance  that 
the  elements  in  the  dictionary  D;  have  in  synthesizing  x;  at  the 
same  sample.  This  method  is  effective,  but  it  can  be  very  time 
consuming.  A  similar  approach  was  also  proposed  in  [28]. 

A  novel  cartoon  and  texture  image  component-based 
restoration  for  SAR  images  was  proposed  in  [29] .  Specifically, 
a  SAR  image,  considered  as  a  function  /,  is  to  be  decomposed 
into  a  sum  of  two  components  f  =  u  +  v,  where  u  represents 
the  cartoon  or  geometric  (i.e.,  piecewise  smooth)  components 
of  /,  and  v  represents  the  oscillatory  or  textured  components. 
The  second  component  essentially  accounts  for  noise  and  the 
texture  elements.  In  the  proposed  technique,  u  is  estimated  and 
is  considered  as  the  restored  (despeckled)  image,  while  the 
textured  components  of  v  are  not  attempted  to  be  recovered. 
This  is  problematic  since  discarding  the  texture  components 
may  result  in  the  loss  of  important  salient  features  in  a  SAR 
image.  Fig.  1  provides  an  example  of  how  an  image  separates 
into  the  structural  and  textured  components  indicating  the 
importance  of  retaining  these  textural  elements. 

Motivated  by  recent  advances  in  sparse  representation- 
based  image  separation  [30],  we  propose  a  similar  separation- 
based  despeckling  method  so  that  the  image  is  decomposed 
into  a  sum  of  piecewise  smooth  and  textured  elements.  Our 
formulation  is  based  on  finding  sparse  representations  of  these 
elements  from  dictionaries  specifically  suited  to  compress 
them,  and  differs  significantly  by  not  just  treating  the  texture 
and  noise  components  as  complements  of  the  cartoon-based 
estimated  image.  By  taking  advantage  of  the  ability  of  sparse 
representations  in  our  scheme  to  estimate,  we  are  able  to  retain 
important  salient  and  textured  features  in  the  final  estimated 
image. 


A.  Organization  of  This  Paper 

Section  II  reviews  the  statistical  models  for  SAR  images. 
In  Section  III,  we  present  the  component  separation  algorithm 
based  on  surrogate  functionals.  In  Section  IV,  we  show  some 
of  the  results  on  both  real  and  simulated  data.  We  present  the 
concluding  remarks  in  Section  V. 


(C)  (d) 


Fig.  1.  Image  separation,  (a)  Original  image,  (b)  Structural  components, 
(c)  Textural  components,  (d)  Structural  +  textural  components.  Note  that  the 
textural  components  in  (c)  contain  some  important  features. 


II.  Statistical  Models  for  SAR  Images 

Let  Y  e  RNxN  be  the  observed  image  intensity,  X  e  RNxN 
be  the  noise  free  image,  and  F  e  RNxN  be  the  speckle  noise. 
Assuming  the  SAR  image  represents  an  average  of  L  looks, 
Y  is  related  to  X  by  the  following  multiplicative  model  [31]: 

Y  =  FX  (1) 

where  F  is  the  normalized  fading  speckle  noise  random 
variable.  It  follows  a  gamma  distribution  with  unit  mean  and 
variance  j-  and  has  the  following  probability  density  function 
(pdf): 

p(F)  =  Y^LLpL~le~LF  (2) 

where  F  >  0,  L  >  1,  and  T(.)  denotes  the  gamma  function. 

The  natural  logarithmic  transformation  converts  the  multi¬ 
plicative  noise  model  (1)  into  an  additive  model 

Y  =  ln(7)  =  ln(F)  +  \n(X) 

=  F  +  X  (3) 

where  F  =  ln(F)  and  X  =  ln(X).  The  pdf  of  the  random 
variable  F  is  given  by  [31],  [32] 

p{F)  =  — L  LLeFLe~LeF.  (4) 

'  7  r(L) 

The  mean  of  F  is  given  by 

E[F]  =  y( 0,  L)-ln(L) 
and  variance  is  given  by 

var(F)  =  y/(l,  L) 


This  article  has  been  accepted  for  inclusion  in  a  future  issue  of  this  journal.  Content  is  final  as  presented,  with  the  exception  of  pagination. 


PATEL  et  al.:  SEPARATED  COMPONENT-BASED  RESTORATION  OF  SPECKLED  SAR  IMAGES 


3 


where  ,  , 

/  d  V+1 

^(M)  =  (  —  J  In  r(0 

is  known  as  the  kth  poly  gamma  function. 

Note  that  there  exist  many  statistical  models  for  SAR 
images  [33],  some  based  on  mathematical  approaches  and 
others  based  on  physical  approaches.  In  this  paper,  we  use 
the  logarithmic  transform  to  convert  the  multiplicative  noise 
to  an  additive  noise  and  take  into  account  the  appropriate 
distribution  model  to  handle  the  noise. 

III.  Image  Component  Estimation 

Let  y,  f ,  and  x  be  lexicographically  ordered  vectors  of  size 
N 2  representing  Y,  F,  and  X ,  respectively.  We  assume  that 
the  SAR  reflectivity  field  is  a  superposition  of  two  different 
signals 

X  =  xp  +  X,  (5) 

where  xp  is  the  piecewise  smooth  component  of  x  and  xt  is 
the  texture  component  of  x.  So  the  additive  model  (3)  can  be 
written  as 


y  =  x  +  f 

=  Xp  +  Xf  +  f .  (6) 

We  further  assume  that  xp  is  compressible  in  a  dictionary 
represented  in  a  matrix  form  as  Dp,  and  similarly,  xt  is 
compressible  in  a  dictionary  represented  in  a  matrix  form 
as  Dr.  Given  Mp,Mt  >  A2,  the  dictionary  Dp  e  Ma^2><mp 
and  Dt  e  RN  xMt  are  chosen  such  that  they  provide  sparse 
representations  of  piecewise  smooth  and  texture  contents, 
respectively.  That  is,  we  assume  that  there  are  coefficient 
vectors  ap  e  Mm^x1  and  at  e  Mm?x1  so  that  xp  =  Dpap 
and  xt  =  D tat.  The  compressibility  assumption  means  that, 
when  the  coefficients  are  ordered  in  magnitude,  they  decay 
rapidly.  The  texture  dictionary  Dt  needs  to  contain  atoms  that 
are  oscillatory  in  nature  such  as  those  found  in  the  discrete 
cosine/sine  transform  and  the  Gabor  transform.  The  dictionary 
Dp  should  be  able  to  process  images  with  geometric  features 
such  as  edges.  The  matrix  D^,  could  represent  some  type  of 
wavelet,  shearlet,  curvelet,  or  contourlet  dictionary. 

One  can  recover  the  SAR  reflectivity  field  x  by  estimating 
the  components  xp  and  xt  via  ap  and  at  by  solving  the 
following  variational  problem: 

ap,at  =  argmin  A||ap||i  +  A\\at\\i  +  yTV(Dpap) 

CLp,(Xt 

+  ^lly-Dp«p -Dr«tll2  (7) 

where  TV  is  the  TV  (i.e.,  sum  of  the  absolute  variations  in 
the  image)  and  for  an  A-dimensional  vector  x,  \\.\\q  denotes 
the  £q- norm,  0  <  q  <  oo,  defined  as 


The  two  components  are  the  corresponding  representations 
of  the  two  parts  and  can  be  obtained  by  xp  =  D pdp  and 
Xf  =  Df<5f.  This  notion  of  separating  a  signal  into  different 


morphologies  using  sparse  representations  is  often  known  as 
morphological  component  analysis  (MCA)  [30]. 

Instead  of  seeking  the  sparse  sets  of  coefficients  ap  and  at 
directly  and  then  inverting  the  representations,  it  is  possible 
to  directly  seek  the  images  whose  transform  coefficients  or 
dictionary  representations  are  sparse.  This  corresponds  to 
solving  the  following  optimization  problem: 

\p,xt  =  arg  mini  ||DNp  111  +  A||Djxf  ||i  +  yTV(xp) 

Xj 9>Xt 

+  f  ||y  -xp  -xt\\22  (9) 

where  D,  denotes  the  Moore-Penrose  pseudo-inverse  of  D,. 
Here,  we  have  assumed  that  the  two  redundant  dictionaries  are 
of  full  rank  and  we  can  obtain  the  analysis  operator  from  the 
synthesis  by  using  the  following  relations: 

ap  =  D \xp 

at  =  Djxr. 

One  of  the  major  advantages  of  using  (9)  is  that  it  requires 
searching  lower  dimensional  vectors  rather  than  longer  dimen¬ 
sional  representation  coefficient  vectors.  This  increases  numer¬ 
ical  efficiency  and  decreases  memory  constraints. 

A.  Iterative  Shrinkage  Algorithm 

Various  methods  can  be  used  to  obtain  the  solution 
of  (7)  [34],  [35].  In  this  section,  we  derive  a  fast  convergent 
iterative  shrinkage  algorithm  by  a  method  of  using  separable 
surrogate  functionals  (SSF)  to  solve  the  separation  problem 
posed  in  (7)  [35]  -  [37].  For  simplicity,  we  assume  that  D  = 
[Dp,Df]  and  discard  the  TV  component  for  the  discussion 
given  here.  The  objective  function  in  (7)  can  then  be  rewritten 

as  i 

/(a)  =  l||a||1  +  -||y  — Da||2  (10) 

where  a  contains  both  the  piecewise  smooth  and  texture  parts. 
Let 

d(a,ao )  =  1 1| a  -  «ol||  -  -||Da  -  Da0|l2  (U) 

where  ao  is  an  arbitrary  vector  of  length  N2  and  the  parameter 
c  is  chosen  such  that  d  is  strictly  convex.  This  constraint  is 
satisfied  by  choosing 

c  >  ||DtD||2  =  lmax(DrD)  (12) 

where  lmax(DrD)  is  the  maximal  eigenvalue  of  the  matrix 

DrD. 

Adding  (11)  to  (10)  gives  the  following  surrogate  function: 

/(a)  =  A||a||i+l||y-Da||2-l-^||a-aoll2-hlD«-D«oll2- 

O3) 

This  surrogate  function  f(a)  can  be  re-expressed  as 

/(a)=A  +  -||a||i  +  d|a— x0||2  (14) 

c  2 

where 

1  T/ 

x0  =  -D  (y  -  Da0)  +  ao 
c 


(15) 
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and  A  is  some  constant.  Let  (a)+  denote  the  function 
max  (a,  0),  and  sign(x)  be  the  signum  function  defined  as 


Given  that 


signO)  = 


-1, 

0, 

1, 


if  a  <  0 
for  a  =  0 
for  a  >  0. 


<5i(x)  =  sign(x)(|x|  -  X)+ 


(16) 


is  the  element-wise  soft-thresholding  operator  with  threshold 
/l,  the  global  minimizer  of  the  surrogate  function  is  given  by 

&sol  =  Sk,c  (xo) 

-Dao)+«o^-  (17) 

It  was  shown  in  [36]  that  the  iterations 

ak+x  =  Sx/c  (;V(y  -  D«b  +  «*)  (18) 

converge  to  the  minimizer  of  the  function  /  in  (10),  where  the 
superscript  k  indicates  that  it  is  the  value  for  the  kth  iteration. 
By  breaking  the  above  iteration  into  the  two  representation 
parts  and  considering  the  TV  term,  we  get  the  following 
iterative  updates  that  essentially  solves  (7): 

«p+1  =  Sx/c  -  DP«P  -  D tdf)  +  dp)  (19) 

ak+1  =  DTpHSyk  (h7^**1)  (20) 

ak+1  =  Sx,c  (Id  J  (y  -  T)pdkp  -  D  tdkt)  +  afj  (21) 

where  H  is  the  undecimated  Haar  wavelet  dictionary. 
A  detailed  description  of  the  undecimated  Haar  wavelet  trans¬ 
form  can  be  found  in  [38].  We  have  replaced  the  TV  correction 
term  by  a  redundant  Haar  wavelet-based  shrinkage  estimate,  as 
this  seems  to  give  the  best  results.  This  adjustment  is  applied 
only  to  the  piecewise  smooth  component  to  control  the  ringing 
artifacts  near  the  edges  caused  by  the  oscillations  of  the  atoms 
in  the  dictionary  Dp.  The  same  adjustment  was  used  in  [30], 
and  the  substitution  was  partially  motivated  by  observing  the 
connection  between  TV  and  the  Haar  wavelet  given  in  [39]. 

The  iterations  presented  above  can  be  extended  to  handle  the 
analysis  formulation  in  (9).  This  is  simply  done  by  modifying 
iterations  (19)— (21)  as  follows: 

xk+1  =  Dp. Sx/c  (^Dp(y  -  Xp  -  xf)  +  DpXp'J  (22) 

xk+1  =  H Syk  (Hrx*+1)  (23) 

xf+ 1  =  Dt  .Sx/c  (V  (y  —  xkp  —  xf )  +  Djx^  .  (24) 

We  summarize  the  algorithm  for  recovering  the  two  sep¬ 
arated  components  of  a  SAR  image  in  Fig.  2.  In  step  3 
of  the  algorithm  in  Fig.  2,  \\.\\oo  denotes  the  ^oo-norm. 
For  an  A-dimensional  vector  x,  it  is  defined  as  Hxlloo  = 
max(|xi|, . . |jcjvD- 


Input:  y,  c. 

Initialization:  Initialize  k  =  1  and  set 

x°  =  0,  x°  =  0,  r°  =  y  —  x°  —  x°,  and  A0  = 

\  (l|Djy||oo  +  l|Dfy||oo)  • 

repeat: 

1 .  Update  the  estimate  of  xp  and  x*  as 

Xp  =  Dp.<SAf=  +  Dtx*-1) 

=  H57,  (HTXp) 

xl\  =  Bt.Sxk  (rfc_1)  +  Djxf-1^)  . 

2.  Update  the  residual  as 

rfe=y-x£— x*. 

3.  Update  the  shrinkage  parameter  as 

Afc  =  I  (HDjrioo  +  ||D^rft|U)  . 

until:  stopping  criterion  is  satisfied. 

Output:  The  two  components  xp  =  x^  and  xt  =  x£. 

Fig.  2.  SSF  iterative  shrinkage  algorithm  to  solve  (9). 


Once  the  two  denoised  components  of  x  are  estimated,  we 
obtain  the  final  estimate  of  x  as 

x  =  exp(x^  +Xf).  (25) 

Since  the  logarithmic  transformation  introduces  a  bias  on 
the  final  estimate,  we  correct  it  along  with  the  exponential 
transformation  [12],  [25]. 


IV.  Experimental  Results 

In  this  section,  we  present  the  results  of  our  proposed 
despeckling  algorithm  and  compare  them  with  the  enhanced 
Lee  filter  [2]  and  some  recent  state-of-the-art  methods 
[25],  [30].  We  also  compare  our  results  with  a  Stein-Block 
thresholding  (SBT)  method  proposed  in  [40].  This  method 
was  shown  to  be  nearly  minimax  over  a  large  class  of 
images  in  the  presence  of  additive  bounded  noise.  This  method 
requires  a  threshold  parameter,  which  we  set  to  the  theoretical 
value  4.505  as  derived  in  [40].  Furthermore,  we  compare  the 
performance  of  our  combined  dictionary-based  approach  to 
despecking  with  that  of  a  fixed  transform-based  despecking 
method.  In  particular,  we  apply  soft  thresholding  on  the 
subband  coefficients  of  the  wavelet  transform.  We  call  the 
resulting  method  wavelet-based  thresholding  (WT).  For  the 
MCA  method  [30],  we  use  the  curvelet  transform  to  represent 
the  piecewise  smooth  component  and  2D-DCT  to  represent 
the  texture  component. 

In  Fig.  3,  we  display  the  test  images  used  for  different 
experiments  in  this  paper.  In  these  experiments,  we  use  the 
relative  error  (RE)  and  the  equivalent  number  of  looks  (ENL) 
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(a) 


(b) 


(c) 


Fig.  3.  Images  used  in  this  paper  for  different  experiments,  (a)  256  x  256 
Cameraman  image,  (b)  512  x  512  Fields  image,  (c)  512  x  512  Nimes  image. 


to  measure  the  performance  of  the  routines  tested,  given  by 

l|x-x||2 
I|X||2 
mean2 
variance 

where  the  mean  and  the  variance  are  measured  within  a 
homogeneous  region  (see  [21]  and  [2]).  A  large  ENL  value 
corresponds  to  better  speckle  reduction. 


RE  = 

ENL  = 


A.  Dictionaries 

In  our  experiments,  we  use  a  dictionary  corresponding 
to  the  shearlet  transform  to  represent  the  piecewise  smooth 
component.  Shearlets  extend  the  traditional  wavelets  by  allow¬ 
ing  waveforms  to  be  defined  not  only  at  various  scales  and 
locations  but  also  at  various  orientations  [41].  The  shearlet 
transform  is  best  suited  for  representing  images  with  edges  and 
anisotropic  structures.  Moreover,  numerical  results  give  evi¬ 
dence  to  the  superior  behavior  of  shearlet-based  decomposition 
algorithms  when  compared  to  curvelet-based  and  contourlet- 
based  algorithms  [41],  [42].  This  is  the  reason  why  we  choose 
shearlets  instead  of  curvelets  and  contourlets  in  our  approach. 
A  brief  discussion  about  the  shearlet  transform  is  given  in  the 
Appendix.  Fig.  4  shows  some  atoms  from  a  shearlet  dictionary. 
In  our  implementation,  we  used  the  nonsubsampled  shearlet 
transform  with  the  decomposition  structure  [3,  3,  4,  4],  which 
determines  the  number  of  directions  in  the  scales  from  coarse 
to  fine.  As  a  result,  the  size  of  is  N2  x  56N2. 

The  discrete  cosine  transform  is  known  to  closely  approx¬ 
imate  the  Karhunen-Loeve  transform  for  a  class  of  random 
signals  known  as  first-order  Markov  processes  that  model 
several  real-world  images.  Explicitly,  the  coefficients  of  the 


Fig.  4.  Atoms  from  a  shearlet  dictionary.  Each  block  represents  the  result 
of  applying  the  shearlet  transform  for  a  particular  scale  and  orientation  after 
applying  it  to  a  centered  impulse  response. 


Fig.  5.  Atoms  from  the  2D-DCT  dictionary.  Each  block  represents  the  result 
of  applying  the  inverse  2D-DCT  transform  to  impulse  responses,  which  are 
centered  at  various  locations. 

DCT  for  an  N  x  N  image  x  =  (xWjW)  are 

2  N-lN-l 

=  *  zz  Xm,n  COS (3t(m  -  1/2)0'  -  l/2)/A0 

m= 0  n= 0 

x  cos(;r {n  -  1/2 ){k  -  1/2 )/N). 

It  is  clear  the  oscillatory  nature  of  the  basis  functions  can 
represent  repetitious  elements  common  to  texture  components, 
which  is  why  they  are  commonly  used  for  this  purpose.  A  few 
atoms  from  the  2D-DCT  dictionary  are  shown  in  Fig.  5.  The 
number  of  atoms  in  the  dictionary  Dr  are  64  x  N2. 

B.  Parameters 

From  the  discussion  in  Section  III,  the  parameter  c  should 
be  chosen  such  that  c  >  /ImaxCD^D^  +  DtDj).  This  can  be 
satisfied  by  choosing  c  >  2.  In  particular,  the  value  we  chose 
was  c  —  3. 

The  Haar  shrinkage  value  yk  in  step  1  of  the  algorithm 
is  3dk  ,  where  ok  is  the  standard  deviation  of  the  noise 

Xp7  Xp 

estimated  by  using  a  median  estimator  on  the  finest  scale  of 
the  Haar  wavelet  coefficients  of  xk  [19]. 
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(c)  (d) 


Fig.  6.  (a)  Noisy  image,  L  =  4,  RE  =  0.498.  (b)  Restored  image  using  our 

method,  RE  =  0.118.  (c)  Noisy  image,  L  =  4,  RE  =  0.500.  (d)  Restored 
image  using  our  method,  RE  =  0.065. 

C.  Stopping  Rule 

We  change  the  threshold  value  of  Xk  during  each  iteration 
according  to 

=l-{\\T?prk\\00  +||Dfrr*||00)  (26) 

and  stop  the  iterations  when  hk  <  To,  where  T  «  2.1  [43]. 
This  way,  we  stop  the  iterations  when  the  residual  is  at  the 
noise  level  and  the  noise  is  rejected  in  each  component.  Note 
that  this  step  of  the  algorithm  assumes  the  noise  variance  to 
be  known.  This  is  not  a  problem  since  the  noise  variance  can 
easily  be  estimated  by  using  a  median  estimator  on  the  finest 
scale  of  the  wavelet  coefficients  of  y  [19]. 

D.  Results  on  Simulated  Data 

In  Table  I,  we  report  the  results  of  experiments  on  the 
simulated  data  with  a  various  number  of  looks.  We  did  not 
have  access  to  the  codes  of  [21],  [25],  and  [22].  Hence,  we 
report  the  RE  values  reported  in  the  corresponding  papers. 
Figs.  6  and  7  show  the  noisy  and  restored  images  for  some 
of  the  experiments  with  the  simulated  data.  It  can  be  seen 
from  these  figures  and  the  results  in  Table  I  that  our  method 
performs  favorably  over  some  of  the  competitive  methods 
for  speckle  reduction.  In  particular,  this  is  the  case  when 
the  number  of  looks  is  greater  than  3.  Furthermore,  these 
results  clearly  indicate  that  an  improvement  is  achieved  when  a 
combined  dictionary  approach  is  used  to  restore  a  SAR  image, 
as  can  be  seen  by  comparing  the  results  of  our  method  with 
that  of  SBT,  MCA,  and  WT  in  Table  I. 

Using  the  Cameraman  image,  in  Figs.  8  and  9,  we  show  the 
evolution  of  the  objective  function  and  RE,  respectively,  as  we 


Fig.  7.  Experiment  with  the  Nimes  image,  (a)  Noisy  image,  L  =  10,  RE  = 
0.315.  (b)  Restored  image  using  our  method,  RE  =  0.163.  (c)  Noisy  image, 
L  =  4,  RE  =  0.501.  (d)  Restored  image  using  our  method,  RE  =  0.207. 
(e)  Noisy  image,  L  =  1,  RE  =  1.001.  (f)  Restored  image  using  our  method, 
RE  =  0.290. 

vary  the  number  of  looks.  Note  that  in  Fig.  9  the  RE  decreases 
significantly  after  five  iterations  and  saturates  around  the  10th 
iteration,  showing  that  the  proposed  method  is  efficient  and 
requires  less  number  of  iterations  compared  to  [21]. 


E.  Results  on  Real  SAR  Data 

In  the  second  set  of  experiments,  we  use  the  real  SAR 
images  shown  in  Fig.  10(a)  and  (b).  These  images  were 
collected  using  the  Sandia  National  Faboratories  Twin  Otter 
SAR  sensor  payload  operating  at  X  band.  Since  the  true 
reflectivity  fields  are  not  available,  we  use  ENF  to  measure  the 
performance  of  our  method.  The  value  of  ENF  is  estimated 
from  the  two  32  x  32  homogeneous  regions  (shown  with  white 
boxes).  We  refer  to  the  region  left  side  of  the  image  as  R1  and 
the  right  side  of  the  image  as  R2.  The  estimated  ENF  values 
are  reported  in  Table  II.  As  can  be  seen  from  this  table,  our 
method  outperforms  the  WT,  MCA,  and  SBT  methods. 
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TABLE  I 

REs  for  Various  Experiments 


Image 

L 

Noisy 

Ours 

SBT  [40] 

MIDAL  [21] 

[25] 

[22] 

WT 

Lee 

MCA 

Cameraman 

13 

0.277 

0.083 

0.104 

0.090 

- 

0.098 

0.098 

0.170 

0.090 

Cameraman 

10 

0.315 

0.090 

0.111 

0.097 

0.091 

- 

0.105 

0.171 

0.101 

Cameraman 

4 

0.498 

0.118 

0.144 

0.124 

0.131 

- 

0.135 

0.178 

0.140 

Cameraman 

3 

0.573 

0.129 

0.156 

0.130 

- 

0.151 

0.147 

0.182 

0.151 

Cameraman 

1 

0.990 

0.184 

0.220 

0.167 

0.192 

- 

0.228 

0.211 

0.210 

Fields 

10 

0.316 

0.054 

0.058 

0.056 

0.055 

- 

0.063 

0.063 

0.052 

Fields 

4 

0.500 

0.065 

0.073 

0.066 

0.066 

- 

0.077 

0.804 

0.071 

Fields 

1 

1.000 

0.099 

0.114 

0.089 

0.096 

- 

0.159 

0.135 

0.110 

Nimes 

10 

0.315 

0.163 

0.312 

0.170 

0.174 

- 

0.195 

0.273 

0.186 

Nimes 

4 

0.501 

0.207 

0.223 

0.217 

0.217 

- 

0.247 

0.277 

0.221 

Nimes 

1 

1.001 

0.290 

0.312 

0.301 

0.314 

- 

0.346 

0.298 

0.320 

TABLE  II 

Estimated  ENL  Values 


Image 

Region 

Original 

WT 

SBT 

Ours 

Lee 

MCA 

Fig.  11(a) 

R1 

3.291 

30.366 

49.028 

79.357 

5.518 

50.23 

Fig.  11(a) 

R2 

3.236 

30.428 

84.936 

133.962 

8.263 

90.89 

Fig.  12(b) 

R1 

4.117 

44.717 

75.041 

146.379 

11.35 

85.75 

Fig.  12(b) 

R2 

3.997 

37.253 

74.641 

101.491 

11.30 

86.70 

x  10' 


Fig.  8.  Objective  function  value  as  a  function  of  the  iteration  number  for 
the  experiments  with  a  Cameraman  image  with  various  numbers  of  looks. 

The  despeckling  results  from  various  methods  are  shown 
in  Figs.  ll(b)-(d)  and  12(b)-(d).  The  ratios  of  the  original 
image  to  the  filtered  images,  referred  to  as  the  noise  images, 
are  shown  in  Figs.  ll(e)-(h)  and  12(e)-(h).  It  is  evident  from 
these  figures  the  single  dictionary-based  reconstructions  such 
as  SBT  and  WT  suffer  from  noticeable  artifacts.  The  MCA 
method  based  on  curvelets  and  2D-DCT  provides  good  recon¬ 
struction,  but  it  removes  a  lot  of  point  targets.  Our  combined 
dictionary  approach  clearly  provides  good  reconstructions  and 
removes  most,  if  not  all,  of  these  artifacts  and  preserves 
point  targets.  Note  that  we  have  not  compared  our  method 
on  the  real  SAR  images  with  the  other  methods  since  their 
implementations  were  not  available  to  us. 

F.  Computational  Efficiency 

In  our  image  separation-based  despeckling  method,  the  most 
computationally  intensive  part  is  in  finding  the  coefficients 


Fig.  9.  RE  versus  number  of  iteration  curves  for  the  experiments  with  a 
Cameraman  image  with  various  numbers  of  looks. 


Fig.  10.  Real  SAR  images  used  for  the  experiments. 

from  a  shearlet  dictionary.  Using  Matlab  on  a  Windows 
system  with  an  Intel  Core  2  CPU  2.16  GHz/3.00  GB  processor, 
one  iteration  of  the  entire  algorithm  takes  around  4.15  s.  On 
average,  our  method  takes  about  45  s  to  process  an  image 
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(e)  (f)  (g)  (h) 

Fig.  11.  Despeckled  images.  Restored  using  (a)  SBT  method,  (b)  MCA  method,  (c)  WT,  and  (d)  our  method,  (e)-(h)  are  the  noise  images  corresponding  to 
the  filtered  images  in  (a)-(d),  respectively. 


(e)  (f)  (g)  (h) 


Fig.  12.  Despeckled  images.  Restored  using  (a)  SBT  method,  (b)  MCA  method,  (c)  WT  and  (d)  our  method,  (e)-(h)  are  the  noise  images  corresponding  to 
the  filtered  images  in  (a)-(d),  respectively. 


of  size  256  x  256.  The  performance  of  our  algorithm  can 
be  enhanced  by  using  a  more  efficient  shearlet  transform 
implementation  which  is  parallelizable.  Each  iteration  of  our 
algorithm  has  the  complexity  of  0(N 2  lo g2(A0). 

G.  Discussion 

Note  that  we  have  applied  our  component  separation  method 
on  the  log-transformed  images.  However,  one  can  directly 
apply  this  method  on  the  SAR  image  without  applying  the  log 


transformation  to  separate  the  piecewise  smooth  component 
and  the  texture  component.  For  instance,  SAR  model  (1)  can 
be  rewritten  as 

y  =  xf 

=  X  +  (f  -  l)x  (27) 

where  l  =  [l,...,l]r.In  this  setting,  the  first  and  the  second 
term  in  (27)  can  be  viewed  as  xp  and  xt ,  respectively.  In  partic¬ 
ular,  if  one  has  designed  dictionaries  specifically  to  compress 
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these  terms,  then  one  can  view  the  first  term  as  the  restored 
image.  However,  designing  dictionaries  to  compress  the  indi¬ 
vidual  terms  in  (27)  is  nontrivial.  One  can  adapt  dictionary 
learning  methods  [44]  to  learn  these  components.  This  in  turn 
would  require  a  collection  of  training  images  containing  only 
the  piecewise  smooth  reflectivity  fields  and  images  containing 
only  signal-dependent  noise.  This  formulation  is  also  relevant 
in  the  case  when  SAR  images  contain  a  strong  additive  noise 
component  where  the  purely  multiplicative  model  may  not  be 
adequate. 

V.  Conclusion 

We  proposed  a  new  method  of  speckle  reduction  in  SAR 
imagery  based  on  separating  an  image  into  various  compo¬ 
nents.  Unique  to  this  approach  is  the  ability  to  use  specific 
dictionaries  of  representations  suited  for  separation  with  an 
iterative  scheme  that  is  able  to  retain  important  features. 
Experiments  showed  that  this  method  performs  favorably 
compared  to  other  competitive  methods.  This  new  process  is 
also  valuable  for  many  SAR  image  understanding  tasks  such  as 
road  detection,  railway  detection,  ship  wake  detection,  texture 
segmentation  for  agricultural  scenes,  and  coastline  detection. 
In  addition,  specific  dictionaries  could  be  designed  to  be  used 
with  this  procedure  to  capture  unique  signatures  while  dealing 
with  speckle  removal. 

Appendix 

Shearlet  Transform 

The  shearlet  construction  can  be  considered  a  natural  exten¬ 
sion  of  wavelets  into  two  dimensions  [45].  Its  representative 
elements  are  defined  by  the  2-D  affine  system 

(V'astto  =  I  det  Mas\A  -  t)  :  t  e  JR2} 

where 

Mas  =  (o  l)  (o  Vo)  (28) 

is  a  product  of  a  shearing  and  anisotropic  dilation  matrix  for 
(a,  s)  G  R+  x  R.  The  generating  function  yj  is  such  that 

to  =  *(&,&)  = 

where  yj\  is  a  continuous  wavelet  for  which  T/i  g  C°°(R ) 
with  supp  Tfi  C  [—2,  1/2]  U  [1/2,  2],  and  y/2  is  chosen  so  that 
^2  e  C°°(M),  supp  ^2  C  [-1, 1],  with  %  >  0  on  (-1, 1), 
and  ||  \ff2 1| 2  =  1-  Under  these  assumptions,  a  function  /  g 
L2(M2)  can  be  represented  as 

f  da 

/(*)=/  /  /  (f,  Vast)  Va&tix)  ^fdsdt 

Jr2  J-oo  7o  a5 

for  a  g  M+,  s  e  M,  and  t  e  M2.  The  operator  SH  defined  by 
SHf(a,s,t )  =  (/,  y/ast) 

is  referred  to  as  the  continuous  shearlet  transform  of  /  g 
L2(M).  It  is  dependent  on  the  scale  variable  a ,  the  shear  s, 
and  the  location  t. 


The  collection  of  discrete  shearlets  is  given  by 
{Wj,e,k  =  I  det  A\i/2  y(Bl  Aj x  -  k)  :  j,l  eZ,kc  Z2}  (29) 
where 

s=(J!)  A=(oJ0-  <30) 

Shearlets  form  a  Parseval  frame  (tight  frame  with  bounds  equal 
to  1)  for  L2(M2),  given  the  appropriate  choice  of  the  generating 
function  yt  (see  [41]  for  details).  An  M- channel  filterbank 
implementation  can  be  done  by  using  the  techniques  given 
in  [46].  As  a  consequence,  its  implementation  has  a  complexity 
of  (D(N2  lo g2(AO)  f°r  an  N  x  N  image. 
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